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ABSTRACT 

New turbulence modeling options recently imple- 
mented for the 3-D version of Proteus, a Reynolds-aver- 
aged compressible Navier-Stokes code, are described. The 
implemented turbulence models include: the Baldwin- 
Lomax algebraic model, the Baldwin-Barth one-equation 
model, the Chien k-s model, and the Launder-Sharma k-e 
model. Features of this turbulence modeling package 
include: well documented and easy to use turbulence mod- 
eling options, uniform integration of turbulence models 
from different classes, automatic initialization of turbu- 
lence variables for calculations using one- or two-equation 
turbulence models, multiple solid boundaries treatment, 
and fully vectorized L-U solver for one- and two-equation 
models. Validation test cases include the incompressible 
and compressible flat plate turbulent boundary layers, tur- 
bulent developing S-duct flow, and glancing shock wave/ 
turbulent boundary layer interaction. Good agreement is 
obtained between the computational results and experi- 
mental data. Sensitivity of the compressible turbulent solu- 
tions with the method of y + computation, the turbulent 
length scale correction, and some compressibility correc- 
tions are examined in detail. The test cases show that the 
highly optimized one- and two-equation turbulence models 
can be used in routine 3-D Navier-Stokes computations 
with no significant increase in CPU time as compared with 
the Baldwin-Lomax algebraic model. 

INTRODUCTION 

Rapid advancements in computer technology and 
numerical algorithms have made it possible to routinely 
perform 3-D Reynolds-averaged Navier-Stokes (RANS) 
calculations for a number of practical turbulent flow prob- 
lems. At NAS A Lewis, a computer code called Proteus has 
been developed to solve the Reynolds-averaged, unsteady 
compressible Navier-Stokes equations. This computer code 
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has been developed with an emphasis on code readability, 
modularity, and documentations. Proteus is available in the 
2-D/axisymmetric 1,2 and 3-D 3 versions. Both versions 
have previously been released with the Baldwin-Lomax 4 
algebraic model (BL) and the Chien 5 k-e model (CH). Tne 
recent implementation effort has added the Baldwin-Barth 6 
one-equation model (BB), the Launder-Sharma 7 k-e model 
(LS), and an automatic starting routine for turbulent calcu- 
lations to the 3-D version of Proteus. In addition, the com- 
pressibility corrections of Zhang et aL 8 , Sarkar 9 , Zeman- 
free shear 9 , Zeman-wall bounded 9 , and Wilcox 9 , as well as 
a turbulent length scale correction, similar to the ones pro- 
posed by Vuong and Coakley 10 and described by Horst- 
man 11 , are available for any of the k-e models above. 

Proteus is a general purpose RANS code that makes 
no prior assumptions on the type or geometry of the fluid 
dynamics problem to be solved. Ideally, it is desirable to 
have a single universal turbulence model in Proteus that 
can handle every turbulent flow problem in that code’s 
domain of applicability. However, none of the present tur- 
bulence models is universal, and the next best thing is to 
have a number of different well proven turbulence models 
available in Proteus as turbulence modeling options. Then, 
the Proteus code users can just try out the various turbu- 
lence models available and select those that work be t for 
their particular applications. 

In this implementation effort, there are four major 
objectives: 

• The turbulence models should complement each 
other in terms of applicability, robustness, and 
accuracy so that together, they can cover the 
widest possible range of turbulent analysis that 
can be done. 

• Although the formulations of the models can vary 
widely between different classes of models, the 
implementation of these models should 
nevertheless be uniform, highly integrated, and 
efficient. 

• The turbulence modeling options need to be 
modular, easy to modify, easy to use, and well 
documented 

• The implemented models should be validated as 
much as possible using standard benchmark test 
cases. 


1 


In this paper, the approaches used in the current work 
to meet the above objectives will be discussed. The new 
Proteus turbulence modeling capability and its features will 
be described. Comparisons will be made between the 
experimental data and the computational results using dif- 
ferent turbulence models for the test cases of incompress- 
ible and compressible flat plate turbulent boundary layers, 
turbulent developing S-duct flow, and glancing shock 
wave/turbulent boundary layer interaction. Sensitivity of 
the compressible turbulent solutions with y + computation, 
turbulent length scale correction, and compressibility cor- 
rections will be examined. Finally, details in using the vari- 
ous turbulence options in 3-D Proteus computations will be 
given. 

TURBULENCE MODELS 

The four turbulence models used in this implementa- 
tion effort can be organized into 3 major classes: algebraic, 
one-equation, and two-equation models. While all of these 
turbulence models are well known, the BL algebraic model 
is perhaps the most popular turbulence model for RANS 
computations. It is well proven, and considerable experi- 
ence in using it has developed over the years. The BB one- 
equation model has recently been developed, and it has 
successfully been applied in a variety of 3-D turbulent cal- 
culations 12 " 14 . The CH and the LS are well proven low- 
Reynolds number k-e models, and they have also been used 
in a number of implicit as well as explicit Navier-Stokes 
codes 15 " 21 . These low-Reynolds number k-e models are 
applicable all the way to the solid wall with no additional 
wall treatments required. 

The LS k-e model is the only turbulence model in 
this group that does not require y + in its formulation. This 
is desirable, since in a generalized 3-D full RANS analysis 
with multiple solid boundaries, it is difficult to compute y + . 

Generally, using more sophisticated turbulence mod- 
els requires more work and turbulence modeling expertise 
from the code user. Using one- and two-equation turbu- 
lence models, for example, will require initial conditions 
and boundary conditions for the turbulence variables in 
addition to those for the mean flow variables. Also, the sta- 
bility requirements of the numerical methods used to solve 
the turbulence equations will need to be considered 
together with the mean flow equations. However, careful 
implementation of these models can significantly reduce 
the amount of work and frustration for the code user, and 
many features in the current turbulence modeling package 
have been implemented for this reason. 

Upwind differencing for the turbulence equations is 
used for both numerical stability and ease of use. Auto- 
matic initialization of turbulence variables is available to 
simplify the task of starting calculations with one- and two- 
equation models. A minimal amount of information is 
required from the user to start a turbulent calculation with 


the Proteus code. The only inputs required for a turbulent 
calculation are the turbulence model type and the boundary 
conditions for the turbulence variables (if a one- or two- 
equation turbulence model was used). User inputs for other 
turbulence parameters or constants are optional. 

Algebraic Model 

The BL algebraic model is the first turbulence model 
implemented into the Proteus code. It is available in both 
the 2-D/axisymmetric and 3-D Proteus versions. Besides its 
primary use as a turbulence model, it is also used as part of 
the automatic initialization procedure for computations 
using one- and two-equation turbulence models. 

Although using the BL algebraic model to compute 
complex 3-D turbulent flows can be questionable, espe- 
cially for separated flows and for flows with multiple shear 
layers and/or multiple solid boundaries, this model is nev- 
ertheless simple, fast, and versatile. Most importantly, it 
gives reasonably good results in all of the test cases consid- 
ered in this report In this model, the formula for the turbu- 
lent viscosity is: 

inner - y crossover) 

^ _ ^ (4) (y >y ) ^ 

outer v J n J crossover-' 

Where y n is the normal distance from the wall and y^ss- 
over is the smallest value of y n at which the inner and the 
outer values of turbulent viscosity are equal. 

The inner turbulent viscosity is defined as 

dinner =P 1 Vl Re r ^ 

with 

1 = Ky n [l-e- y:/A ’] (3) 


| co| is the magnitude of the vorticity and 



where 

Pi can either be p w or p L 
V] can either be v w orv L 

The subscripts w and L in the above terms denote 
wall and local conditions, respectively. The reason for the 
factor Re r in eq. (2) and eq. (4) and in the equations that 
follow in this report is that all of the variables have been 
nondimensionalized in Proteus. x w is computed as the 
product of the fluid molecular viscosity at the wall and 
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either the normal gradient of the tangential velocity at the 
wall or | ©| at the wall. |co| is slightly cheaper to compute 
than the normal gradient of the tangential velocity, so it is 
used in most of the computations. For 2-D boundary layer 
flows, y„ + computed using | co| is the same as using the nor- 
mal gradient of the tangential velocity. For other flows con- 
sidered in this report, the error in y n + computed from | ©| is 
small and did not significantly affect the final solutions. 

Options are available in Proteus to use the Launder- 
Priddin length scale modification for eq. (3) and to use the 
Spalding and Kleinstein inner turbulent viscosity formula 
instead of eq. (2), but none of these was used in the present 
investigation. 

The outer turbulent viscosity is defined as 

(ff t ) outer = KC C pPFw akc F Kll ,bR e r (5) 


A + = 26.0; Cep =1.6; Ciaeb = 0.3; 

C wlc = 1.0; K = 0.0168; k = 0.4 

In Proteus, modifications to the BL model were made 
so that it can also be used for turbulent free shear flows. 
These modifications are given in detail in the Proteus man- 
ual sets 1 ' 3 , and they will not be described here, since all of 
the test cases considered in this investigation involve only 
wall bounded flows. 

One-Equation Model 

The BB one-equation model is a recently developed 
model for the turbulence quantity k^/e that avoids the need 
for an algebraic turbulent length scale. Also, this model has 
less numerical difficulties than the k-e models. Experiences 
with the BB model in Proteus calculations showed that it is 
very similar to the BL model in terms of speed, robustness, 
and grid point clustering requirement 


with 


The actual implementation of this model uses the fol- 
lowing nondimensionalized equations: 


F = v F 

wake J max max 

F max in eq. (6) is the maximum value of 
F(y n ) = y|co[(l-e- y:/A ’) 

and yma* is the value of y n corresponding to F max . Note 
that eq. (6) and eq. (7) only apply to wall bounded flows. 

As recommended by Degani and Schiff 22 , F max is 
taken to be the first peak of F(y„) when searched from the 
solid wall. To prevent the selection of spurious peaks of 
F(y„) near the wall, a peak is considered to have been 
found when the value of F(y„) drops below FPMIN*F max , 
where FI MIN is a user adjustable factor. The default value 
for FP1.ICN in Proteus is 0.9. 


3w 3F dG 3H 
•=r- +^-+xr- + xr- = S 
d t dx dy d z 

(9) 

where 


w = pa 

(10) 

1 da 
F = pua- 5 j-n. s 

(11) 

G = pva- 5— 

p Re, ‘ l dy 

(12) 

1 3a 

H = pwa- - — u. 
p Re ^*dz 

(13) 


The function Fj^eb in eq. (5) is the Klebanoff inter- G : 

mittency factor. The modified FjQgb formula used in Pro- S = (c^fj ~ c ei) 
teus is given by: 1 


^Kleb C*. m j n + 




1 +B 


^KIeby n 

\ ^max 


-1 


( 8 ) 


Qc,min ^ a user adjustable constant normally. When 
this constant is set to zero, eq. (8) becomes the usual 
Klebanoff expression. However, when the BL model is 
used to generate the initial turbulent viscosity field to start 
the k-e calculations, C kmin is set to be around 0.1. This will 
prevent the turbulent viscosity values from becoming too 
small in the uniform freestream regions and help minimize 
starting problems with the k-e computations. The constants 
used are: 



k 2 

a = — 
e 

1 _ , ^ 
- ( c e2 c el) ~~T 
°£ K 

C7 r = 0.25 CT C 

p t = CppaDjDj 


(14) 

(15) 

(16) 

(17) 

(18) 

(19) 
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c e2 V c e2 J Ky 


( 20 ) 

( 21 ) 

(22) 


transport equations for the turbulent kinetic energy k and 
the isotropic dissipation rate e. Other k-e models use the 
transport equation for the total dissipation rate e, 23-25 , but 
the wall boundary condition for the isotropic dissipation 
rate e is simpler (e^, = 0), and numerical experiments dur- 
ing this investigation have found that the CH and the LS 
models are more forgiving numerically than the k-E models 
using e l , especially in the near wall region. 

The LS k-e model does not require y n + , and this is an 
important advantage for complex 3-D flows with multiple 
solid boundaries. On the other hand, the near wall formula- 
tion of the CH model is simpler and more robust numeri- 
cally. The nondimensionalized k-e equation in vector 
notation is the same as eq. (9), where 


(23) 

(24) 
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F = 
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Lpe_ 
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(26) 


G = 
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1 3e 


, 1 3kn 
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1 3e 
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(30) 


(31) 


F - 
1 Re r 

|G)| 

3 s 


< 27 > S = 

K = 0.41; 

c E l = 1.2; 

Cg2 = 2.0; 


Cji = 0.09; 

A + = 26; 

& 

+ 

II 

; O 


Bj = 0.4; 

B 2 = 0.01; 

e 0 = 10' . 

^ = 


P k — Re r pe + D + C k 

e £ 2 

Ci p k - - Re r c 2 f2P^- + E + C E 


Mt 

a. 


(32) 


(33) 


In eq. (15), the term with the subscript (ij,k) is held 
constant at a point (i ode) and not differenced like the rest of 
Pa. That term is actually outside of the partial derivatives, 
but it is included inside Pa for programming convenience. 
The turbulence variable a = k 2 /e is nondimensionalized by 
v„ a reference kinematic viscosity value, and the other 
mean flow variables are nondimensionalized using stan- 
dard reference values as described in the Proteus manuals. 






= C^p- 


(34) 


(35) 


Two-Equation Models 

The CH and LS models are two-equation low-Rey- 
nolds number k-e models. They are well proven, and a lot 
of experiences have been gained on the use of these two 
model over the years. Both the CH and LS models use the 


~ rT P ‘ Pi 3 pk?2 


(36) 


where Pj and P 2 are given by eq. (23) and eq. (24), respec- 
tively. 
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For the CH model, the following damping functions and 
constants apply: 


f H = 1 - e~° 011Sy * 

(37) 

f 2 = 1 - ?e _Rc,1/36 

(38) 


where 


c a = 1.35; c 2 =1.80; c tl = 0.09; 

<7 k =1.00; ct £ =1.30 

For the LS model, the following damping functions and 
constants apply: 


f 2 = 1 - 0.3e~ Re ' 1 (41) 

Where Re l is defined in eq. (39), and 

Cj = 1.44; c 2 = 1.92; c R = 0.09 
a k = 1.00; a e = 1.30 

D and E in eq. (32) are the extra near wall modifications 
used by the CH and the LS models, and they are given by: 


f 2pk 

Re r y 2 „ 


(CH) 

~ (^-Vk) 

l Re r dXj 

dXi 

(LS) 

r 2pe -y/2 


(CH) 

Re r y l C 



f A ) 

(LS) 

Re 2 p WXjdxJ 

{ dxjdxJ 


The term D is present because the isotropic dissipa- 
tion rate e is used. According to Jones and Launder 26 , term 
E is present to improve the computed peak level of turbu- 
lence kinetic energy at y + = 20. The turbulence variables k 
and e have been nondimensionalized by u r 2 and p r u r 4 /ji. r , 
respectively. 

C k and C e in eq. (32) are the compressibility correc- 
tions, and the definitions for those terms will depend on the 
particular compressibility correction used. 


Five compressibility correction options for the k-e 
models were examined: (1) a general correction due to 
Zhang et al. 8 (ZH), (2) the Sarkar correction (SA), (3) the 
Zeman correction for free shear flow (ZF), (4) the Zeman 
correction for wall bounded flow (ZW), and (5) the Wilcox 
correction (WI). Options number 2-5 were implemented as 
described by Wilcox 9 . 

The Zhang et al. generalized compressibility correc- 
tion was obtained by including the Sarkar correction and a 
number of other models for the exact compressible terms in 
the k-e model equations. For this implementation, these 
terms are: 

C k = -P k a 1 M t 2 -Re r pe t (a-a 2 )M t 2 


1 ( ~\ dp 9P 

Re r [p 2 CT p J 3 x i ax i 


1 f K yp a<T ij 
Re 2 Ap 2 cr p J ax i ax j 

(44) 

c e = pe t P 2 [- | - n CY- !) + |c,] 

(45) 

where 


s 

•* N) 

II 

* £ 
H ^ 

(46) 

( du ■ 3uA o 
cr.. = p + ^-2- -fpP,6.. 

» l3x ; 3x- 3^ 2 » 

(47) 


a = 0.5; aj = 0.4; a 2 = 0.2; 

n = 0.7; CT p = 0.5 

The remaining four compressible corrections only 
apply to the k equation so that C e is zero for these correc- 
tions. For these compressibility corrections, C k is defined 
as follows 

C k = -Re r pe t ^F(M t ) (49) 

Sarkar’s Model (SA) 

K = 1 (50) 

F(M l )=M [ 2 (51) 

Zeman’s Model (ZF & ZW) 
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4* = 0.75 


F(M t ) = 


1 - exp - 


(M t -M,o) 


2.1 


(52) e because there is more uncertainties in the model- 
ing of the e variable. In general, the constants C TL and 
y max + will be different for different flows, and the determi- 
nation of these constants will require some knowledge 
about the turbulent length scale in the near wall region of 
the flow under consideration. 


0 

s 

1 

X 


(53) 

where H(x) is the Heaviside step function, and 


Mto = 0.10; A = 0.60 

for free shear (ZF) 


Mjo = 0.25; A = 0.66 

for wall bounded (ZW) 

Wilcox’s Model (WI) 



4* = 1-5 


(54) 

F(M t ) = [Mf-MfjHfM,- 

MJ 

(55) 

Where Mj 0 is 0.25 




According to Coakley and Huang 27 , the dissipation 
rate e in standard k-e models can collapse abnormally near 
reattachment points for separated flow and high speed 
flows. This leads to an erroneous strong growth in turbulent 
length scale, resulting in unrealistically high skin friction 
and heat transfer rates. To correct this defect, a turbulent 
length scale correction similar to the ones proposed by 
Vuong and Coakley 10 and described by Horstman 11 is also 
available for both of the k-e models. With this correction, 
the turbulent length scale in the wall region is limited so 
that it can never be greater than a constant times the dis- 
tance from the wall (C-j^ x y n ). The default value for the 
constant C-j^ is 2.5, which corresponds to the von Karman 
constant of k = 0.4. This length scale correction is only 
applied to a region near the wall. In this work, a maximum 
value for y n + is used to define the upper limit of this near 
wall region. A good discussion in the limit for the near wall 
region was given by Rodi 28 . The default maximum y n + 
value used for the wall region is y max + = 20. The correction 
is applied as follows: 


Y + Computation for 3-D geometries 

The determination of the distance from a solid wall is 
difficult for complex 3-D wall bounded flows with arbitrary 
geometries. In Proteus, the calculation of y n and y n + is 
done in a stand-alone, fully vectorized subroutine. The user 
can either manually indicate the solid portions of the com- 
putational boundaries (point-by-point, if necessary) or let 
Proteus automatically identifies them from the no-slip wall 
boundary condition. After the solid portions of the compu- 
tational boundaries have been identified, the straight line 
distances from all of the solid portions of the computa- 
tional boundaries to an interior point are computed. Then 
the distances are compared, and a nearest solid wall is 
selected. Finally, y n and y n + are computed with respect to 
that nearest solid wall. If there are no solid wall present, 
then y n and y n + are set to a very high numerical value. This 
has the effect of driving all of the Van Driest style turbulent 
damping terms with the form [1 - exp(-y n + /A + )] to 1, which 
is expected for turbulent flows away from solid walls. 
Since the subroutine is separate from the rest of the pro- 
gram, the calculation of y n and y n + can be readily modified 
to accommodate any custom geometry. 

Applying the above algorithm to a hypothetical 2-D 
grid in fig. 1 , one can see that for an interior point A, the y n 
value at A is simply the straight line distance between A 
and B, and the y n + value at A will be computed using the 
wall properties at B. For point C, y n is the straight line dis- 
tance between C and D. 

Numerical Algorithm for One- and Two-Equation 
Models 

The partial differential equation represented by eq. (9) 
is a scalar equation for the BB model and a vector equation 
for the k-e models. The solution procedures for the BB 
model and the k-e models are analogous. Therefore, only 
the k-e solution procedure is discussed below. 


e = 


max e. 


,3/2 


Re r (C TL 



(56) 


Applying the generalized grid transformation to eq. 
(9) and putting it into the conservative form results in 


for y + < y + 

J n J max 

With this lower limit imposed on e, the turbulent 
length scale will not be larger than the quantity C^l x y n in 
the specified near wall region. Horstman 11 limited the tur- 
bulent length scale by imposing an upper limit on the tur- 
bulent kinetic energy. However, in this study, the turbulent 
length scale is limited by limiting the collapse of the dissi- 
pation rate e instead. This was done because of the above 
observation of Coakley and Huang 27 regarding the collapse 


3w 3F 3G 3H * 
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H = He — H d - H M1 — H M2 

(60) 
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( 68 ) 
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(72) 


„ The diffusion terms F M1 , F M2 , G^ , Gm 2 , H m1 , and 
H M2 have mixed second-order derivatives, and they are 
lagged in time on the right hand side to keep the block 
matrix tridiagonal. 

An L-U approximate factorization numerical scheme 
described by Hoffmann 29 was used. The approximately 
factored equations for a two sweep L-U scheme are 


Upward sweep 

e.Ax 
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0, Ax 
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= A w* 


( 73 ) 


(74) 


Where. A, R, C,J), E, P, and Q are, the Jacobian matri- 
ces of F c , F d , G^, Gq, He, H d , and S , computed as 
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where 
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(82) 

(83) 

(84) 

(85) 

(86) 

(87) 

( 88 ) 

(89) 


JC e 2 = P 2 [-|-n(Y-l)+|c 1 ] (97) 

And for the other compressibility corrections 


JC kl = JD£*F(M,) 

(98) 

JC k2 = -Re£*F(M t ) 

(99) 

O 

(n 

II 

O 

(100) 

*-t 

O 

N> 

II 

O 

(101) 


In eq. (73) and eq. (74), the Jacobian matrices A, C, 
and E are differenced using the first-order upwind differ- 
encing, and the Jacobian matrices B, D, and P are differ- 
enced using the second-order central differencing. In the 
RHS, the first-order upwind differencing is used to approx- 
imate the convective terms and the second-order central 
differencing is used for the diffusive terms. The central dif- 
ferencing operators for the Jacobian matrices B, D, and P 
have been split into a forward differencing part and a back- 
ward differencing part so that, for example: 


a^B = O^B)' + 0 5 B) + (102) 

(90) 

where (3-B) is the backward differencing part of 
the central differencing operator 3,B , and (5^B) + is the 
forward differencing part. This is done so that only back- 

(91) ward differencing operators are present in the upward 
sweep, eq. (73), and only forward differencing operators 
are present in the downward sweep, eq. (74). 


In the solution algorithm, the upward sweep is fol- 
(92) lowed by the downward sweep. When either eq. (73) or eq. 
(74) is applied at each interior point of the computational 
domain, only one unknown per equation needs to be 
solved. Therefore, eq. (73) and eq. (74) are solved point by 
point 


(93) Assuming that the computational block has a dimen- 

sion of (nl,n2,n3). The right hand side of eq. (73) is known 
at the current time level n. To advance the solution to the 
time level n+1, the upward sweep is first marched starting 
from the lower comer of the computational block, point 
(2,2,2) to the upper comer of the computational block, 
point (nl-l,n2-l,n3-l). During the upward sweep, eq: (73) 
is solved for the intermediate unknown Aw* at a point 
( ' (il,i2,i3) using data at points (il-1 ,i2,i3), (il,i2-l,i3), and 

(il,i2,i3-l). This is possible because the left hand side of 
eq. (73) contains only backward differencing operators. 
(95) Then the downward sweep is marched in the opposite 

direction of the upward sweep, from point (nl-l,n2-l,n3-l) 
to point (2,2,2). During the downward sweep, eq. (74) is 
(9g) solved for the final unknown Aw at a point (il,i2,i3) 
using data at points (il+l,i2,i3), (il,i2+l,i3), and 
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(il,i243+l). This is possible because the left hand side of 
eq. (74) contains only forward differencing operators. 
Finally the solution is advanced to the new time level using 

w = w + A w (103) 

The marching order in both sweeps can be manipu- 
lated to achieve better efficiency for the L-U solver. For 
example, a straight forward Fortran do-loop for the upward 
sweep is as follows: 


do 

10 

il = 2,nl-l 

do 

10 

i2 = 2, n2-l 

do 

10 

i3 = 2,n3-l 


★ 

★ 

10 continue 

The above do-loop will give the correct results, but it 
does not take full advantage of the vectoiization capability 
of a vector computer. Since the loop has two nested do- 
loops, only the innermost loop is vectorized. On the other 
hand, if the marching is done in the direction normal to the 
diagonal planes il+i2+i3 = constant, then a do-loop can be 
constructed with only one nested do-loop as follows: 

do 10 iplane = l,nplane 
do 10 ipoint = 1, npoint 
★ 

★ 

★ 

10 continue 

Where nplane is the number of diagonal planes in the 
3-D computational block, and npoint is the number of grid 
points contained a diagonal plane. In general, npoint will 
vary from plane to plane. It turns out that the points (il- 
I,i2,i3), (il42-l,i3), and (il,i2,i3-l) are all located in the 
plane (iplane-1), and they are known (for an upward 
sweep). As the result, the only nested inner loop in the 
above do-loop can be vectorized over every point in a diag- 
onal plane. 

To implement the above marching scheme, an 
addressing scheme is needed to translate the do-loop vari- 
ables (ipoinffiplane) to the grid indices (il,i2,i3) so that the 
flow properties at (il4243) can be conveniently retrieved in 
the inner do-loop during the computation of the block 
matrices. There are many ways that this can be accom- 
plished, and, in this project, a simple scheme has been 
devised to compute the il, i2, and i3 grid indices from the 
do-loop variables ipoint and iplane. This scheme does not 
require any machine-dependent routines, and will work for 
any Fortran 77 compiler. Basically, this scheme works as 
follows: 

1. Outside of the inner loop, the il index of every 
point in the current diagonal plane iplane is stored 


in the array iloc(ipoint). 

2. Outside of the inner loop, the diagonal line index 
of every point in a diagonal plane is computed and 
stored in the array line(ipoint). 

3. Inside the inner loop, the il , i2, and i3 grid indices 
can be computed from the iloc(ipoint) and 
line(ipoint) arrays as follows: 

11 = Hoc (ipoint) 

12 = -il + line (ipoint) + 3 

13 = iplane - il - i2 + 5 

Arrays iloc(npoint) and line(npoint) are necessary for 
a fully vectorized inner loop. This addressing scheme is 
illustrated in fig. 2. In this diagram, a computational block 
of (5x6x4) is used. The points on the computational bound- 
aries are assumed to be known through the application of 
the explicit boundary conditions, and the unknown interior 
points are indicated by dots of different shapes. The dots 
with the same shapes belong to the same diagonal plane. It 
can be seen that for any interior point (i 14243) belonging 
to a diagonal plane iplane, the points (il- 14243), (il42- 
l,i3), and (il42,i3-l) belong to the diagonal plane iplane-1, 
and the points (il+l,i2,i3), (il42+143), and (il,i2,L3+l) 
belong to the diagonal plane iplane+1. This makes it possi- 
ble to vectorize the above L-U marching scheme. 

With some effort, the code can be fully vectorized in 
the inner do-loop using the above marching scheme, result- 
ing in a very fast L-U solver for the one- and two-equation 
models. This scheme can also be used in a Navier-Stokes 
solver. 

Boundary values for k and e are computed explicitly 
after the downward sweep. If spatially periodic boundary 
condition in the ^ direction is specified, then an additional 
set of grid points is added at il = nl+1, and flow properties 
at il = nl+1 are explicitly set to the values at il = 2. With 
periodic boundary condition in the il direction, the upward 
sweep is marched from point (2,2,2) to point (nl,n2-l,n3- 
1), and the downward sweep is march from point (nl,n2- 
l,n3-l) to point (2,2 ,2). An analogous procedure is used for 
periodic boundary conditions in the other two spatial direc- 
tions. 

Automatic Initialization Procedure for One- and 
Two-Equation Models 

With the BL algebraic model, the user only needs to 
specify the initial conditions for the Navier-Stokes equa- 
tion. However, with one- and two-equation turbulence 
models, the user also needs to specify the initial conditions 
for the turbulence variables to start the time marching pro- 
cess. An automatic initialization procedure for the turbu- 
lence variables would greatly simplify the task of the code 
user. This automatic initialization procedure must provide a 
unified, transparent integration of turbulence models, even 
if they are from different classes. Using this automatic 
starting procedure, the user should be able to use any of the 
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turbulence models in. a new or restart calculation without 
worrying about the initial conditions for the turbulence 
variables. 

In Proteus, the type of the turbulence model used in a 
calculation is saved together with a minimum number of 
turbulence variables (zero for algebraic models, one for 
one-equation models, and two for two-equation models) in 
a restart file. For the BL model, no turbulence variable is 
saved. For the BB model, k^yfe is saved, and for the k-e 
models, k and £ are saved. The automatic initialization pro- 
cedure will depend on the types of turbulence model 
selected for the current and previous calculations. 

BL model 

If the BL model is used in the current calculation, 
then jq is calculated using eq. (1). No further initialization 
is needed. 

BB model 

For a new calculation or a restart calculation that uses 
the result of a previous laminar calculation, the turbulent 
viscosity (q is first calculated using the BL model. Then the 
initial value for k^/e is calculated using eq. (19). 

For a restart calculation that uses the result of a previ- 
ous turbulent calculation, the turbulent viscosity |q is first 
calculated using the saved turbulence variable(s) and the 
appropriate expression for |q. Then the initial value for k 2 /e 
is calculated using eq. (19). If the previous calculation used 
the BB model, then only |q needs to be computed. 

CH or LS k-e model 

For a new calculation or a restart calculation that uses 
the result of a previous laminar calculation, the turbulent 
viscosity iq is first calculated using the BL model. Then 
using the assumption of local equilibrium and the CH 
damping function f^, initial profiles for k and £ are gener- 
ated with 



For a restart calculation that uses the result of a previ- 
ous turbulent calculation, the turbulent viscosity |q is first 
calculated using the saved turbulence variable(s) and the 
appropriate expression for |q. If the previous calculation 
used either the BL or BB models, then eq. (104) and eq. 
(105) are used to obtain the initial profiles for k and £. 
However, if either the CH or the LS model was used in the 
previous calculation, then the k and £ profiles saved from 
the previous calculation are used as the initial profiles. 

The previous automated procedure has been found to 
perform well for the test cases considered. For example. 


fig. 3 shows that the automatically generated profiles for k 
are actually very close to the converged k profiles for the 
cases of flat plate turbulent boundary layer, turbulent devel- 
oping S-duct flow, and glancing shock wave/turbulent 
boundary layer interaction. For the flat plate case, the near 
wall peak in the initial k profile matches the experimental 
data even better than the converged k profile. The local 
equilibrium assumption, eq. (104), and eq. (105), can also 
be used to generate the inflow profiles for k and e when 
such profiles are not otherwise available. 

It is recommended that the turbulent calculations 
using one- or two-equation turbulence models be started 
from a converged turbulent solution. This will at least pro- 
duce some reasonable initial profiles for the turbulence 
variables, minimizing the starting problems. In addition, 
the initial CFL number used for starting a k-£ calculation 
should be well below one. After reasonably smooth profiles 
for k and £ have been achieved, then the CFL number can 
be increased. 

VALIDATION TEST CASES 

Validation test cases include the incompressible and 
compressible flat plate turbulent boundary layers, turbulent 
developing S-duct flow, and glancing shock wave/turbulent 
boundary layer interaction. These test cases represent a 
good mixture of compressible and incompressible wall- 
bounded flows. The grid sizes for these cases range from 
20,655 grid points for the incompressible flat plate case to 
321,776 grid points for the glancing shock wave/turbulent 
boundary layer interaction case. All computations were 
done on a Cray YMP, and the CPU time required is 
reported for each test cases. In all of the test cases below, 
the default artificial viscosity options for the Navier-Stokes 
equations in Proteus were used. Since all of the flows con- 
sidered are steady, convergence was assumed to be reached 
when either the residuals or interested flow properties 
stopped changing with more iterations at the highest CFL 
number used in a calculation. 

All of the input and output parameters used in Proteus 
must be nondimensionalized by the appropriate reference 
conditions, and it is important to specify the correct refer- 
ence conditions in a Proteus calculation. For this reason, all 
relevant reference conditions for a particular calculation 
are discussed below. 

Incompressible Flat Plate Turbulent Boundary 
Layer 

Problem Description 

Incompressible turbulent flow over a flat plate at zero 
pressure gradient with a freestream Mach number of 0.3 
was computed using the BL, BB, CH, and LS turbulence 
models. The computational domain for the flat plate ranges 
from Re x = l.OxlO 4 to Re x = 10xl0 6 . Comparisons are 
made with the Klebanoff 30 turbulent boundary layer pro- 
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files at Re 0 = 7700 and the Wieghart 31 skin friction data. 
Reference Conditions 

For convenience, the reference Reynolds number Re r 
was chosen to be IxlO 6 . To simulate the incompressible 
flow, the freestream Mach number of 0.3 was used. The ref- 
erence Mach number was set equal to the freestream Mach 
number. Assuming standard sea level values of air tempera- 
ture, kinematic viscosity, and pressure for the correspond- 
ing reference values, the reference length was calculated to 
be 0.47763 lft 


Grid Geometry 

This is a 2-D problem, but a 3-D 81x51x5 grid (fig. 
4a) was used with 5 stations in the third dimension. The x 
(streamwise) coordinate for the flat plate under consider- 
ation runs from x = 0.01L r to x=10.0L r with 81 evenly dis- 
tributed grid points. This streamwise distribution of grid 
points allows the x station corresponding to the Klebanoff 
profiles to be in the middle of the computational domain 
and away from the inflow and outflow boundaries, mini- 
mizing the effects of boundary conditions on the solution at 
that station. 51 grid points are used in the y (normal) direc- 
tion. The lower y coordinates at the plate surface (lower 
boundary) are simply y = 0 for all x. For the upper 
freestream boundary, the maximum y coordinate increases 
with the streamwise distance along the plate to keep 
approximately the same number of grid points within the 
boundary layer as it grows. This is analogous to a boundary 
layer code that adds more grid points in the normal direc- 
tion as the boundary layer grows. To obtain the maximum y 
coordinates, the turbulent boundary layer thickness at an x 
station was estimated from the power law formula 32 : 


8 = 


0.38x 


Re 


1/5 


(106) 


In the flat plate calculations, approximately 90% of 
the grid points are placed inside the turbulent boundary 
layer in the normal direction, and the rest of the grid points 
are placed in the freestream region. To resolve the near wall 
region, grid points are packed near the plate surface using 
an algebraic formula so that the y + value of the nearest 
point to the wall is approximately from 0.1 at the leading 
edge to 0.9 at the trailing edge. At the x station correspond- 
ing to the Klebanoff experiment, the nearest y + is about 0.5, 
and there are approximately 10 points within y + < 10.0. 

Initial Conditions 

To start the BL calculation, the initial velocities u and 
v were calculated using the Blasius solution for the laminar 
flat plate boundary layer. The non-dimensionalized pres- 
sure, P/Pp was set to 1.0 everywhere. Calculations using 
the BB and CH models were done using the BL solution as 
the initial conditions. The initial values of k 2 /e, k, and e for 
the BB and the CH model were obtained using the auto- 
matic starting procedure described above. The LS calcula- 
tion was started from the CH solution. 


Boundary Conditions 

For all of the calculations for this test case, constant 
stagnation enthalpy was assumed. Therefore, no boundary 
condition was needed for temperature. 

At the plate leading edge, u and v were specified 
using the Blasius solution. This was done because the Re* 
value at that location is well within the laminar region, w 
was set to zero. The pressure was set to the freestream 
value. E, gradients for k 2 /e, k, and e were set to zero for the 
first approximately 800 iterations. After that, k 2 /e, k, and e 
profiles were held fixed for the rest of the calculations. This 
was done to smooth out any kinks or sharp angles in the 
upstream profiles of k^/e, k, and e that were produced when 
the BL model was used in the automatic initialization pro- 
cess. 

At the plate trailing edge, u, v, k, and e were linearly 
extrapolated in the ^ direction, w and the E, gradient of k 2 /e 
were set to zero. The pressure was set to the freestream 
value. 

At the flat plate surface, the velocity, k, e, k^/s, and 
the normal pressure gradient were set to zero. 

At the outer boundary, u and P were set to the 
freestream values. The velocity component w and the T] 
gradients of v, k, and e were set to zero. The nondimen- 
sional value of k 2 /e was set to 0.5. 

Finally, the £ gradients of all flow variables were set 
to zero to simulate 2-D flow. 

Computational Details and Results 

For the BL and the BB calculations, the starting CFL 
number was 5, and it was gradually increased to a final 
value of 20 during the calculations. For the CH calcula- 
tions, the starting and ending CFL values were 1 and 20, 
and for the LS calculation, the minimum and maximum 
CFL numbers were 1 and 15. LS calculations with CFL 
number large, than 15 diverged. 

For the BL calculation, the skin friction distribution 
stopped changing after 3839 iterations. The BB, CH, and 
LS calculations took 6631, 7531, and 8632 iterations, 
respectively, to reach the same level of convergence. The 
CPU times required for the BL, BB, CH, and LS calcula- 
tions are 6.031xl0' 5 , 6.650xl0' 5 , 6.822xl0' 5 , and 
8.667xl0' 5 sec/iteration/grid point, respectively. As can be 
seen, the differences in the CPU times required for compu- 
tations with different turbulence models are relatively 
small, even with models from different classes. For exam- 
ple, the CH model only required 13% more CPU sec/itera- 
tion/grid point than the BL model for this test case. Note 
that the required CPU times above are for full Navier- 
Stokes calculations with no simplifying assumptions for 
the boundary layer flow. 

The plots of the velocity profile, fig. 5a and fig. 5b, 
show that all of the turbulence models accurately compute 
the turbulent velocity profile as compared with the Kleban- 


11 


off data at Ree = 7700. The differences between the models 
are quite small. In fig. 5c, the Reynolds shearing stress pre- 
dictions by all of the turbulence models are also in very 
good agreement with the Klebanoff profile with the BB and 
BL models giving better agreement for the Reynolds shear- 
ing stress in the outer region of the boundary layer, and the 
CH and LS models giving better prediction of the peak of 
the Reynolds shearing stress in the near wall region. In fig. 
5d, the CH model is seen to predict the peak in k better than 
the LS model, but in the outer region, both models give 
practically the same k profile. The k 2 /e profiles obtained 
using the BB, CH, and LS models are compared in fig. 5e. 
There were differences between the BB model and the k-e 
models, but the general trend and magnitude are approxi- 
mately the same for all models. 

In fig. 5f, the local skin friction coefficient predictions 
by the CH and LS are slightly better than the BL and BB 
models at higher values of Ree, but the differences are rela- 
tively small between the different turbulence models. 

The LS is well known for not accurately predicting 
the peak in k in the near wall region. However, in this par- 
ticular calculation, that defect does not seem to affect its 
accuracy in velocity, Reynolds shearing stress, k 2 /e, and 
skin friction predictions as compared with the CH model. 
The LS and the CH results are almost identical for those 
mean flow properties. 

Compressible Flat Plate Turbulent Boundary 
Layer 

Problem Description 

Compressible turbulent flow over a near adiabatic flat 
plate at zero pressure gradient with a freestream Mach 
number of 3.0 was computed using the BL, BB, CH, and 
LS turbulence models. The computational domain for the 
flat plate ranges from Re x = 3.7xl0 4 to Re x = 2.2xl0 7 . 
Comparisons are made with the Van Driest II correlation 
(VDH) as described by Hopkins and Inouye 33 . In addition, 
velocity and temperature profiles at the x station corre- 
sponding to Re 0 = 9751 are also compared with the Kim 
and Settles 37 inflow profiles for their glancing shock wave/ 
turbulent boundary layer interaction experiment. 

Reference Conditions 

The solution of this calculation is used as the inflow 
boundary condition for the glancing shock wave/turbulent 
boundary layer interaction calculation. Therefore, the flow 
conditions for this calculation have been carefully chosen 
to match those from the Kim and Settles experiment 34 ' 37 . 
For convenience, the reference length L r was set equal to 
the incoming boundary layer thickness at Re 0 = 975 1 in the 
Kim and Settles experiment, which is 0.009908 ft. The ref- 
erence Mach number was set equal to the freestream Mach 
number of 3.0. Using the isentropic relation, the freestream 
static temperature was calculated to be 189 deg. R, and this 
value was used as the reference temperature. Assuming the 


reference pressure to be equal to the freestream static pres- 
sure, the reference density was calculated to be 0.04667 
lbm/ft 3 . The reference Reynolds number was calculated 
from the above reference values to be 1.8694 x 10 5 . 

Grid Geometry 

This is a 2-D problem, but a 3-D 91x71x5 grid (fig. 
4a) was used with 5 stations in the third dimension. The x 
(streamwise) coordinate for the flat plate under consider- 
ation runs from x = 0.2L,. to x=120.0L,. with 91 evenly dis- 
tributed grid points. This streamwise distribution of grid 
points allows the x station corresponding to the Kim and 
Settles inflow profiles to be approximately in the middle of 
the computational domain and away from the boundaries, 
minimizing the effects of boundary conditions on the solu- 
tion at that station. 7 1 grid points are used in die y (normal) 
direction. Following the procedure outlined above for the 
incompressible flat plate calculation, the height of the grid 
is increased in the streamwise direction to keep approxi- 
mately the same number of grid points in the boundary 
layer. To resolve the near wall region, grid points are 
packed near the plate surface so that the y + value of the 
nearest point to the wall is about from 0.02 at the leading 
edge to 0.6 at the trailing edge. At the x station correspond- 
ing to tne Kim and Settles inflow profiles, the nearest y + is 
about 0.4, and there are approximately 14 points within y + 
< 10 . 0 . 

Initial Conditions 

To start the BL calculation, the initial velocities u and 
v were calculated using the Blasius solution for the laminar 
flat plate boundary layer, w was set to zero, and P/P r was 
set to 1.0 everywhere. The boundary layer temperature pro- 
file was calculated using the following relation 38 : 


2 

T = T + (T -T ) — - — 
* 1 " w 'u 2c 

c p 

(107) 

where 


rU* 

T - T + 

* w e 2(C p ) e 

(108) 

T w = 1.06T aw 

(109) 

r = 0.896 

(HO) 


Note that eq. (109) was taken from the Kim and Set- 
tles data. 

Calculations using the BB and CH models were done 
using the BL solution as the initial conditions. The initial 
values of k 2 /e, k, and e for the BB and the CH models were 
obtained using the automatic starting procedure described 
above. The LS calculation was started from the CH solu- 
tion. 
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Boundary Conditions 

For all of the calculations in this test case, the energy 
equation was solve simultaneously with the momentum 
and the continuity equations. The molecular viscosity and 
the thermal conductivity coefficients were calculated as 
function of the local temperature using Sutherland’s for- 
mula for air 38 . 

At the flat plate leading edge, u and v velocity compo- 
nents were specified using the Blasius solution, w was set 
to zero. The pressure was set to the freestream value. 
Inflow static temperature profile was held fixed from the 
initial conditions. 1; gradients for k 2 /e, k, and e were set to 
zero for the first approximately 300 iterations. After that, 
k^/e, k, and £ profiles were held fixed for the rest of the cal- 
culations. This was done for the same reason as discussed 
above in the incompressible flat plate test case. 

At the plate trailing edge, £, gradients of all flow vari- 
ables were set to zero. 

At the flat plate surface, the velocity, k, e, k 2 /e, and 
the normal pressure gradient were set to zero. T was set to 
the wall temperature reported in the Kim and Settles exper- 
iment 

At the freestream boundary, u, P, and T were set to the 
corresponding freestream values. The velocity component 
w and the r| gradients of v, k, and e were set to zero. The 
nondimensional value of k 2 /e was set to 0.5. 

Finally, the £ gradients of all flow variables were set 
to zero to simulate 2-D flow. 

Computational Details and Results 

The minimum CFL numbers used are 2.0 for the BL 
model and 0.5 for all the other models. The maximum CFL 
numbers used are 10 for the BB model and 5 for all the 
other models. Note that the higher CFL value of 10 might 
also work for models other than the BB model, but this has 
not been tried due to time constraints. 

For the BL calculation, the skin friction distribution 
remained constant after 2855 iterations. The BB, CH, and 
LS calculations took 2745, 3781, and 3174 iterations, 
respectively, to reach the same level of convergence. The 
CPU times required for the BL, BB, CH, and LS calcula- 
tions are 9.042xl0' 5 , 9.436xl0" 5 , 9.836xl0' 5 , and 
11.77xl0' 5 sec/iteration/grid point, respectively. More 
CPU time is required for this calculation than the incom- 
pressible flat plate calculation, because the energy equation 
was solved. Even so, there are still very small differences 
in the CPU times required for different turbulence models. 
In this test case. The CH model only took 9% more CPU 
lime than the BL model. 

The skin friction plot, fig. 6a, compares the BL, BB, 
CH, and LS skin friction results with the VDII. Also shown 
is a bounding bar denoting a ±10% deviation from the 
VDII. It can be seen that all of the models give results that 
are within this bounding bar. The BL, BB, and the CH 


models produce the best agreement, whereas the LS mod- 
el’s prediction is about 10% lower than the VDII. It should 
be noted that with the BL, BB, and the CH models, the way 
that y + is computed affects the skin friction predictions, and 
this effect will be examined later. In these calculations, y + 
was computed using eq. (4) with Vj = v L and = p w . 

In fig. 6b, the momentum thickness predictions of the 
BL, BB, CH, and LS models are compared with the VDII. 
Again, the LS model results deviates the most from the 
VDII. Note that if y + for the CH model was computed with 
V[= v w and pj = p^ the CH skin friction and momentum 
thickness results will be very similar to the LS results. The 
LS results are not sensitive to y + computation, because that 
model does not use y + in its formulation. 

The velocity and static temperature profiles are com- 
pared with the Kim and Settles data at Reg = 9751 in fig. 6c 
and d. The experimental data were collected as the incom- 
ing boundary layer profiles for a glancing shock wave/tur- 
bulent boundary layer interaction experiment All of the 
turbulence models used give results that are in excellent 
agreement with the experimental data. Also shown in the 
velocity plot is the Schlichting 39 empirical correlation of 
the velocity profile in a compressible boundary layer for a 
freestream Mach number of 2.4. 

Computing y + for incompressible flow is straight for- 
ward. But for compressible flow, there are some ambigu- 
ities in the v and p values which are needed to compute y + 
(see eq. (4)). Baldwin and Lomax 4 recommended Vj=v w 
and Pi=p w in their original report. But Vj=v L was used in 
more recent works in compressible turbulence modeling 6, 
8 . Furthermore, Coakley and Huang 27 reported that the 
choice in pj has important implications for the predictions 
of skin friction and heat transfer, especially for cold walls. 
Therefore, to examine the sensitivity of the BL, BB, and 
the CH results with Vj and p t choices in the y + computa- 
tion, flat plate computations were made using: 

1. The BL model with (Vj= v L , pj=p w ), (v!=v L , 
Pi=Pl). (vi=v w , Pi=p w ), and (v^v^p^Pl) 

2. The BB model with (V!=v L , Pi=p w ) and (V]=v w , 
Pl=Pw) 

3. The CH model with (Vj=v L ,Pi=p w ) and 
(v I =v w ,p I =p w ). 

For the BL model, fig. 7a and b show that the choice 
in V] affects the results more significantly than p ; , and that 
the best results were obtained using (vj= v L , Pi-p w ). 
Using (vpv^ Pi=p w ) as recommended by Baldwin and 
Lomax gives flat plate results that are more that 10% lower 
than the VDII. 

Since the BL results show that pj does not signifi- 
cantly affect the results for these calculations and that the 
best results were obtained using pi=p w , Pi=p w was used 
for the BB and CH models, and only the choice in V] was 
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examined for these models. 

It can be seen in fig. 7c and d that the BB is not as 
sensitive to the choice in Vj as the other turbulence models. 
Both choices in v j produce excellent agreement with the 
VDn. However, using Vi=v w gives slightly lower skin fric- 
tion and momentum thickness results than Vi=v L . 

On the other hand, as shown in fig. 7e and f, using 
Vj=v w in the CH model significantly lowers the skin fric- 
tion and momentum thickness results. In fact, the CH 
model with Vj=v w and pp=p w gives results that are almost 
identical to the LS model. This might be due to the facts 
that the CH and LS models are very close to each other in 
the actual formulations, and that both were derived and cal- 
ibrated in the incompressible turbulent flow regime. 

The unmodified k-e models have been found to fail to 
predict the observed decrease in the spreading rate for the 
compressible mixing layer, and compressibility options are 
available to correct this deficiency. However, these com- 
pressibility corrections can adversely affect the flat plate 
results. As the result, Wilcox 9 recommended using the 
Zeman compressibility option with M to and A set to 0.25 
and 0.66, respectively to minimize its impact on the flat 
plate computation. Wilcox also suggested his own com- 
pressibility correction. To study the effects of these com- 
pressibility corrections on the Proteus flat plate results, the 
Mach 3.0 flat plate computations with the CH model were 
done with and without these corrections. Effects of these 
corrections on the results with the LS model are expected 
to be similar to the CH model, since both models are very 
similar in form. As can be seen in fig. 8a and b, the WI and 
ZW corrections have no effect on the flat plate results, 
while the ZH, SA, and ZF corrections significantly lower 
the skin friction and the momentum thickness Reynolds 
number predictions. 

Turbulent Developing S-Duct Flow 

Problem Description 

Turbulent incompressible flow in an S-duct was com- 
puted using first the BL model and then the CH k-e model. 
The CH calculation was started from a BL converged solu- 
tion. The S-duct geometry and comparison data were 
obtained from an experiment conducted by Taylor et al. 40 
The S-duct in that experiment consists of two 22.5 degree 
bends with constant area square cross section. Turbulent 
experimental data for the S-duct is available for the flow 
Reynolds number of 40,000 (based on the bulk velocity 
and the duct hydraulic diameter, D H ). 

Reference Conditions 

The calculations were done at a flow Mach number of 
0.2 (based on bulk velocity) in order to minimize com- 
pressibility effects and, at the same time, achieve reason- 
able convergence rate with the Proteus code. For 
convenience, the reference Mach number was also set at 


0.2, and the reference Reynolds number was chosen to be 
40,000. Assuming standard sea level values of temperature, 
kinematic viscosity, and static pressure for the reference 
values of TpV,, and P r respectively, the reference velocity 
and the reference length were calculated to be 223.32 ft/sec 
and 0.028658 ft, respectively: 

Grid Geometry 

The computational grid for the S-duct is shown in fig. 
4b. The grid was created using the Gridgen grid generation 
software package for the Iris graphics workstation. The 
computational grid ranges from 7.5 D H upstream of the 
start of the S-bend to 7.5 Dh downstream of the end of the 
S-bend with 81x31x61 grid points in the x, y, and z direc- 
tions, respectively. Since the S-duct is symmetric with 
respect to the y = constant plane, only half of the duct is 
discretized. To resolve the viscous layer, grid points were 
tightly packed near the solid walls so that the y + value of 
the first grid point off of the solid walls is a little less than 
0.5. 

Initial Conditions 

To start the BL calculation, the initial static pressure 
was set equal to the reference pressure, and the fluid veloc- 
ity components u,v, and w were set to zero everywhere in 
the duct. 

To start the CH calculation, the initial values of u, v, 
w, P, and m for the S-duct flow were obtained from the con- 
verged BL solution. The initial profiles for k and e were 
obtained using the automatic initialization procedure 
described above. 

Boundary Conditions 

For both the BL and the CH calculations, constant 
stagnation enthalpy was assume, eliminating the need for 
solving the energy equation and specifying its initial and 
boundary conditions. 

At the duct inlet, the total pressure was specified to be 
Po/P r = 1.028281121. This was computed using the isen- 
tropic formula with the freestream static pressure and the 
bulk Mach number. Zero streamwise gradient was specified 
for u velocity component. The v and w velocity compo- 
nents were set to zero. Zero streamwise gradients were 
specified for k and e for the first 41 iterations of CH calcu- 
lation. After that, the inlet k and e profiles were held fixed 
for the rest of the computation. Around 20-50 iterations at 
the lowest CFL number used is usually enough to smooth 
out those profiles without significantly changing the magni- 
tudes of either k or e. 

At the duct exit, the duct exit pressure was specified 
to be Pe/P r = 0.98415512. This was found by trial and error 
in order to match the experimental mass flow rate in the 
duct. Zero streamwise gradient was specified for all veloc- 
ity components, k, and e. 

At the solid walls, fluid velocity, k, e, and normal 
pressure gradient were set to zero. Standard symmetry 
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boundary condition was used in the symmetry plane. 


Computational Details and Results 

The BL calculation was done with the following 
sequence of CFL numbers 


CFL 

1.0 

5.0 

10.0 


No. of iterations 
100 
200 

After 300 iterations 


And for the CH calculation. 


CFL 

1.0 

3.0 

10.0 


No. of iterations 

762 

1445 

After 2207 iterations 


For the BL calculation, the mass flow rate throughout 
the duct was found to be constant to within 0.1% after 6278 
iterations. Started from the BL solution, the CH calculation 
took another 3667 iterations to achieve the same uniformity 
in the mass flow rate. The CPU times required for the BL 
calculation and the CH calculation are 5.168xl0' 5 and 
6.506xl0' 5 sec/iteration/grid point, respectively. In this cal- 
culation, the CH model required approximately 26% more 
CPU time than the BL model. 

The S-duct streamwise velocity contours and the 
coordinate system used in plotting the results are shown in 
fig. 9. The coordinate system r-z used here is the same as 
that defined in the original experiment. In fig. 10, the wall 
pressure predictions of the BL and CH calculations are 
compared with the experimental data. Note that the refer- 
ence pressure used to calculate c p in the plot was the cen- 
terline wall pressure at streamwise distance = -1.0 D H (one 
hydraulic diameter upstream of the start of the S-bend), r = 
0.5, and z = 1.0. The experimental data shows a larger scat- 
tering of the wall pressures at the streamwise station = -1.0 
Dh than the Proteus result Otherwise, the agreement is 
very good. Both turbulence models correctly computed the 
trend as well as the pressure drop along the S-duct walls. In 
fig. 11 and fig. 12, plots of experimental and computational 
streamwise velocity profiles on the symmetry plane (z = 0) 
and the midspan plane (r = 0.5) of the S-duct are shown for 
the 5 streamwise stations along the S-duct. The agreement 
with the experimental data is generally good for both turbu- 
lence models. The asymmetry in the symmetry plane veloc- 
ity profiles due to the pressure induced secondary motion 
are correctly predicted by the Proteus code. No other turbu- 
lence model was used in the S-duct calculation, because the 
results obtained so far are fairly insensitive to the turbu- 
lence models used. 


Glancing Shock Wave/Turbulent Boundary Layer 
Interaction 

Problem Description 

Glancing shock wave/turbulent boundary layer inter- 
action produced by Mach 3.0 flow past an 10 degree sharp 
fin mounted on a flat plate was computed using the BL and 
CH models. Experimental data are available from the Kiu 
and Settles experiment 34 ' 37 at the Penn State Gas Dynam- 
ics Laboratory. The data available for the purpose of com- 
parison include measured inflow turbulent boundary layer 
velocity and temperature profiles, flat plate surface pres- 
sure, skin friction, and surface flow angles. The surface 
skin friction was experimentally measured by laser inter- 
ferometry. 

Reference Conditions 

This computation uses the same reference conditions 
as described in the compressible flat plate turbulent bound- 
ary layer test case above. 

Grid Geometry 

The geometry for the computational grid used is 
shown in fig. 4c. The grid was generated using the Gridgen 
grid generation software package for the Iris graphics 
workstation. A total of 52x68x91 grid points were used in 
the x, y, and z directions, respectively, a smaller number of 
grid points was used in the x direction than the other two 
directions, because numerical experimentation has found 
that the results are relatively insensitive to grid point distri- 
bution in the x direction, as long as Ax is not too large com- 
pared with the incoming boundary layer thickness. The 
computational domain extended from 6 L r upstream of the 
fin leading edge to 45 downstream of the fin leading 
edge. Recall that L,. was set equal to the incoming turbulent 
boundary layer thickness at Re e = 9751 as reported in ref. 
37. The width of the upstream and downs tre: m boundaries 
are 4.2 L r and 41.3 L n respectively. The upstream width 
was selected so that the computational boundaries are out- 
side of the upstream influence region, and the downstream 
width was selected so that the inviscid shock wave would 
pass through the downstream computational boundary. The 
height of the computational domain is 7 L,, Uniform grid 
spacing of was used in the x direction with Ax = 1.0 L r 
Grid points were clustered in both the y and z directions in 
order to resolve the viscous layers on the flat plate and on 
the fin solid surfaces. The minimum y + value in the y direc- 
tion (on the flat plate) is approximately 0.025 to 0.05, and 
the minimum y + value in the z direction (on the fin) is 
approximately 0.1 to 0.9. Maximum grid spacings in the y 
and z directions are 0.616 L r and 0.531 L p respectively. Six 
streamwise grid points were placed upstream of the sharp 
fin. 

Initial Conditions 

To start the BL calculation, the converged BL fiat 
plate turbulent boundary layer profiles at the inflow x sta- 
tion of the flow velocity, pressure, and temperature were 
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used everywhere except near the fin surface. Near the fin 
surface, a 1/7 th power law velocity profile and a linear tem- 
perature profile were used to smoothly blend the fin surface 
no-slip boundary condition to the rest of the initial flow 
field. 

The converged BL solution was used to start the CH 
calculation. The initial profiles for k and 8 were obtained 
using the automatic initialization procedure as described 
above. 

Boundary Conditions 

At the upstream boundary, the static pressure was set 
equal to the freestream static pressure. The static tempera- 
ture T as well as the velocity components u, and v specified 
at the upstream boundary were interpolated from the con- 
veiged solution of the BL flat plate calculation at the same 
x distance downstream of the Kim and Settles inflow pro- 
files. The w velocity component was set to zero. Zero 
stream wise gradients were specified for k and e for the first 
20 iterations of the CH calculation. After that, the inlet k 
and e profiles were held fixed for the rest of the CH compu- 
tation. 

At the downstream boundary and the far field bound- 
aries (boundaries that are opposite from the fin and the flat 
plate), zero gradients was specified for all flow variables. 

At the flat plate and fin solid surfaces, fluid velocity, 
k, e, and normal pressure gradient were set to zero. The 
static temperature was set to 1.06 T aw as reported in ref. 37. 

Finally, symmetry boundary condition was set at the 
symmetry plane upstream of the fin leading edge. 

Computational Details and Results 

The BL calculation was done with the following 
sequence of CFL numbers 


CFL 

No. of iterations 

0.5 

425 

1.0 

425 

2.0 

850 

3.0 

850 

5.0 

After 2550 iterations 

And for the CH calculation. 

CFL 

No. of iterations 

0.2 

700 

1.0 

700 

3.0 

700 

4.0 

700 

5.0 

700 

Note that the use of CFL values higher than 5.0 was 


possible for both the BL and CH models, but at CFL num- 
bers greater than 5.0 and with the default artificial viscosity 
parameters in Proteus, flow properties started to oscillate 
around shock wave in the inviscid region, even though the 
surface flow properties would still be smooth. Convergence 


was assumed to be reached when the wall pressure, skin 
friction, and surface flow angle distributions no longer 
change with more iterations at the highest CFL number. 
With the above schedule for CFL numbers, the BL calcula- 
tion took approximately 4500 iterations to converge. 
Started from the BL model, the CH calculation required 
approximately another 3500 iterations to converge. The 
CPU times required for the BL calculation and the CH cal- 
culation were 7.888xl0' 5 and 9.757xl0‘ 5 sec/iteration/grid 
point, respectively. The CH k-e model took about 24% 
more sec/iteration/grid point than the BL algebraic model 
for this test case. 

Starting from the converged CH solution, the turbu- 
lent length scale correction took 1353 iterations to con- 
verge, and the Sarkar compressibility correction took 2848 
iterations. Both of these corrections required only slightly 
more CPU sec/iteration/grid point than the standard CH 
computations. Note that when any of the turbulence 
options are changed in a calculation, such as a change in 
the turbulence model types, length scale corrections, or 
compressibility corrections, it is a good idea to reduce the 
CFL number to less than 1.0 to avoid starting problems. 
The CFL number can then be increased gradually after- 
wards. 

When these computations were first attempted, both 
of the BL and the CH computations would always diverge 
in the vicinity of the symmetry plane in front of the fin 
leading edge, regardless of the time step size or the artifi- 
cial viscosity options used. A careful examination of the 
Proteus computations revealed that the extremely tight 
clustering of grid points normal to the symmetry plane has 
caused Proteus to compute unreasonably high gradients of 
the mean flow velocity normal to the symmetry plane, 
causing both the total vorticity magnitude and the produc- 
tion of the turbulent kinetic energy to be too high. This 
caused the turbulent viscosity coefficient in that region io 
be unreasonably high, preventing the BL and CH calcula 
tions from convergence. 

The grid clustering normal to the symmetry plane is 
not good, because there is no real physical reason for it. 
However, it has to be there, because of the structure of the 
grid used. Therefore, to alleviate this problem and allow 
the computations to proceed, the total vorticity magnitude 
and the production rate for the turbulent kinetic energy 
upstream of the fin were computed only in the x-y plane. 
This modification is reasonable, since upstream of the fin, 
the flow is essentially a supersonic 2-D turbulent boundary 
layer flow. For the rest of the flow field, those turbulent 
properties were computed using the full 3-D expressions. 

The BL computed surface streamlines and the general 
layout of the test case is shown in fig. 13. The upstream 
influence and the primary separation lines can clearly be 
seen. From comparisons of the surface flow angles with the 
experimental data below, these computed streamlines are 
quite accurate. 

The BL results are shown in fig. 14. It can be seen that 
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that the skin friction result, fig. 14b, is sensitive to the 
method of y + computation, and the best agreement was 
obtained with the computation of y + using Vj = v r The 
surface pressure and flow angles distributions are seen to 
be relatively insensitive to the method of y + computation. 
Note that the BL Cf prediction shows a small initial rise in 
Cf at P = 30°. Through numerical experimentation, it was 
found that if the number of grid points in the z direction is 
not sufficient, then the BL calculation will not pick up this 
flow feature. On the other hand, the CH k-e calculation will 
always predict this initial rise in c f independent of the num- 
ber of grid points in the z direction. Overall, the agreement 
between the BL results and the experimental data is very 
good. 

The CH results and the effects of the turbulent length 
scale correction are shown in fig. 15. Note that the CH cal- 
culations were done with the method of y + computation 
that gave the best results for the BL calculation above. The 
agreements in the surface pressure and flow angle distribu- 
tions are quite good and, for the most part, are relatively 
insensitive to the turbulent length scale correction applied. 
The standard CH k-e model does quite well in predicting 
the general shape and trend of the Cf curve. In particular, 
the initial rise in c f is always predicted by the CH model 
regardless of the number of grid points used in the z direc- 
tion. This is not the case with the BL model above. How- 
ever, the CH model significantly overpredicts the skin 
friction in the interaction region at b = 15°. This can be 
caused by either the collapse in the turbulent dissipation 
rate as discussed in ref. 27, or insufficient grid-clustering in 
the inner layer of the turbulent boundary layer as reported 
in ref. 41. 

Awa et al. 41 have found that the overshoot in the skin 
friction prediction by the CH model in the reattachment 
region can be corrected by refining the grid in the near wall 
region. D : fferent number of grid points as well as different 
grid clusterings in the y direction have been tried, but they 
did not correct the overshoot in skin friction prediction in 
any significant way. It should be noted that the total number 
of grid points in this calculation was limited by practical 
considerations, and it is possible that the skin friction pre- 
diction of the CH model will improve with further grid 
refinement 

On the other hand, the turbulent length scale correc- 
tion and the compressibility corrections will directly affect 
the skin friction results without requiring more grid points 
or more CPU time. Both of these corrections have the net 
effect of lowering the skin friction prediction. As can be 
seen in fig. 15, the turbulent length scale correction signifi- 
cantly lowers the skin friction as well as increasing the sur- 
face flow angles in the interaction region. A number of 
different values for the turbulent length scale constant C-n_ 
and y + max were tried, but the CH skin friction predictions 
with the length scale correction are still too low (see fig. 
15b) compared with the experimental data. Picking the cor- 
rect turbulent length scale constant and the near wall region 
requires knowledge about the turbulent length scale of the 


flow and necessitates a trial and error approach if this data 
is not available. It is possible that good results can be 
achieved with the right combination of constants. 

Wilcox 9 has found that the Sarkar and other com- 
pressibility correction will improve the k-e skin friction 
prediction in the reattachment region of separations 
induced by shock wave/turbulent boundary layer interac- 
tions. To investigate the effects of the compressibility cor- 
rections in the present problem, the Sarkar correction was 
used. From fig. 16b, it can be seen that the standard Sarkar 
correction brings the CH skin friction prediction in the 
interaction region down to match the experimental data 
precisely. As can be seen in fig. 16a and c, the effects of the 
Sarkar correction on the surface pressure and flow angle 
distribution are small. 

In fig. 17, the effects of the turbulent length scale cor- 
rection and the Sarkar compressibility correction on the 
turbulence variables k, e, and the turbulent length scale 
constant C-^ at a point in the interaction region are shown. 
It can be seen that both of these corrections lower the val- 
ues of k, e, as well as the turbulent length scale constant 

In fig. 17c, for example, the Sarkar compressibility 
correction smoothly lowers the entire Cjl curve, whereas 
the length scale correction will clip the turbulent length 
scale constant Cjl values to below 3.0 for y + values less 
than 20, as specified. It is interesting to note that the k and e 
profiles (fig. 17a-b) computed with the length scale correc- 
tion are still smooth, even though the turbulent length scale 
constant itself is clipped quite abruptly by the length scale 
correction. Also, the turbulent length scale constant tends 
to a value around 2 - 2.5 at y+ > 100. This corresponds to a 
von Karman constant of approximately 0.4 - 0.5. 

Finally, the quasiconical property of this glancing 
shock wave/turbulent boundary layer interaction flow as 
reported by Kim and Settles 42 is investigated. Since Pro- 
teus provided the solution for every point in the computa- 
tional domain, it is easy to compute the skin friction at any 
location on the flat plate. In fig. 18, the best BL and CH 
results of surface pressure, skin friction, and surface flow 
angle distributions are compared for three different radii 
from the fin leading edge. These three radii have been cho- 
sen so that they will be away from both the initial non-con- 
ical inception zone and the downstream computational 
boundary. Note that these radii should really be measured 
from the conical virtual origin which is slightly upstream 
of the fin leading edge, as described in detail by Settles and 
Kimmel 43 . But in this study, the radii were measured from 
-the fin leading edge to facilitate comparisons with the 
experimental data. Moreover, the error in locating the coni- 
cal origin is small for larger values of the radius R. As can 
be seen, both of the BL and CH results show only small 
variations with different radii, and the surface pressure, 
local skin friction coefficient, and surface flow angle distri- 
butions are essentially independent of the radius R. 
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CONCLUSION 

A general turbulence modeling capability has been 
implemented in Proteus, a 3-D compressible Reynolds- 
averaged Navier-Stokes computer code. This new turbu- 
lence modeling capability consists of four turbulence mod- 
els: the Baldwin-Lomax algebraic model, the Baldwin- 
Barth one-equation model, the Chien k-e model, and Laun- 
der-Sharma k-e model. 

Features of the current turbulence modeling imple- 
mentation include: 

• Well documented, modular, and easy to use 
turbulence modeling options. 

• Uniform integration of turbulence models from 
different classes. 

• Multiple solid boundaries treatment. 

• Automatic initialization of turbulence variables 
for calculations using one- and two-equation 
turbulence model. 

• Fully vectorized L-U solver for one- and two- 
equation models. 

Validation test cases showed good agreement between 
the Proteus calculations using different turbulence models 
and the experimental data. The compressible turbulent 
solutions have been found to be sensitive to the method of 
y + computation, turbulent length scale correction, and 
compressibility corrections. 

Finally, the test cases demonstrated that the highly 
optimized one- and two-equation turbulence models can be 
used in routine 3-D Navier-Stokes computations with no 
significant increase in CPU time as compared with the 
Baldwin-Lomax algebraic model. As an example, the 
Chien k-e model only required 9%-26% more CPU sec/ 
iteration/grid point than the Baldwin-Lomax algebraic 
model, depending on the test case. 
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Fig. 1 Computation of y n for multiple solid boundaries 
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Example: nl - 5, n2 - 6. n3 - 4 computational grid 


nplane - nl + n2 + n3 - 8 -7 

ndiag - nl + n2- 5 - 6 

For a point at ipoint - 4, (plane - 5, 

11 - iloc(ipoint) - 3 

12 - -iloc(ipoint) + line(ipoint) +3 - 4 

13 - i plane +5-11-12 -3 


Fig. 2 Grid addressing scheme for the implementation of the L-U marching 
algorithm 
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Fig. 3 Automatic initialization of the k-e models 
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Fig. 5 Incompressible turbulent boundary layer results 
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a Local skin friction distribution 
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Fig. 6 Compressible flat plate turbulent boundary layer results 
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Fig. 7 Effects of method of y n + computation on compressible flat plate turbu- 
lent boundary layer results 
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a. Skin friction, CH model 



b. Momentum thickness, CH model 


Fig. 8 Effects of k-e compressibility corrections on compressible flat plate 
turbulent boundary layer results 
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Flow Direction 


Fig. 9 Streamwise velocity contours and the general geometry for the turbu 
lent developing S-duct test case 
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Fig. 1 0 Wall pressure distribution, turbulent developing S-duct flow 
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Fig. 11 Streamwise velocity distribution on the symmetry plane 
— turbulent developing S-duct flow 
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Fig. 12 
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Streamwise velocity distribution on the midspan plane 
— turbulent developing S-duct flow 
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Fig. 13 Computed surface flow streamlines on the flat plate and the general 
geometry — Glancing shock wave / turbulent boundary layer interac- 
tion 


33 


Surface streamline angle, + (deg.) SurfaCG waH P ressure ' P„ ' 1 P. 






Fig. 14 Glancing shock wave/turbulent boundary layer interaction results — 
The BL turbulence model with effects of y + computations 
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c. Surface streamline angle, CH model 


Fig. 15 Glancing shock wave/turbulent boundary layer interaction results — 
The CH k-e model with effects of the turbulent length scale correction 
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Fig. 1 6 Glancing shock wave/turbulent boundary layer interaction results — 
The CH k-e model and the effects of the Sarkar compressibility cor- 
rection 
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Fig. 17 Glancing shock wave/turbulent boundary layer interaction results — 
Effects of the Sarkar compressibility and length scale corrections on 
some turbulence quantities 
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Fig. 18 Glancing shock wave / turbulent boundary layer interaction results — Quasiconical similarities in computed flow 
properties 
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